FULL AND HALF GILBERT TESSELLATIONS WITH RECTANGULAR 

CELLS 



o ' * 

a 

(N 

5? 
CM 



JAMES BURRIDGE, RICHARD COWAN, ISAAC MA 



Department of Mathematics. University of Portsmouth, Portsmouth, UK. james.burridge@gmail.com 
**School of Mathematics and Statistics, University of Sydney, NSW, 2006, Australia, rcowan@usyd.edu.au 
***Lee Wai Lee Vocational Studies Institute, Hong Kong, isaacma@vtu.edu. hk 



Abstract. We investigate the ray-length distributions for two different rectangular versions of 
Gilbert's tessellation [3]. In the full rectangular version, lines extend either horizontally (with 
east- and west-growing rays) or vertically (north- and south-growing rays) from seed points 
which form a Poisson point process, each ray stopping when another ray is met. In the half 



^ C*| rectangular version, east and south growing rays do not interact with west and north rays. 



For the half rectangular tessellation we compute analytically, via recursion, a series expansion 
for the ray-length distribution, whilst for the full rectangular version we develop an accurate 
simulation technique, based in part on the stopping-set theory of Zuyev [5], to accomplish the 
same. We demonstrate the remarkable fact that plots of the two distributions appear to be 
identical when the intensity of seeds in the half model is twice that in the full model. Our 
' paper explores this coincidence mindful of the fact that, for one model, our results are from 

a simulation (with inherent sampling error). We go on to develop further analytic theory for 
00 ■ the half-Gilbert model using stopping-set ideas once again, with some novel features. Using 

C****" 1 our theory, we obtain exact expressions for the first and second moment of ray length in the 

, half-Gilbert model. For all practical purposes, these results can be applied to the full-Gilbert 

model — as much better approximations than those provided by Mackissack and Miles [4]. 

o 

(N 

1. Introduction 

Consider a stationary Poisson point process in the plane, of intensity A. The particles of this 
^ 1 process are called seeds, aptly so because at a given time t = they each initiate the growth of a 

line. The directions of the lines are randomly distributed, uniformly on (0,7r], and independent 
of each other and of the seed locations. Each line grows bidirectionally from its seed at the same 
rate; thus two rays grow from each seed. When a ray encounters a line that has already grown 
across its path, the growth of that ray stops. Eventually all rays stop growth and a tessellation 
of the plane is formed. 

The completed structure has become known as the Gilbert tessellation after Edgar N. Gilbert. 
It is notoriously difficult to analyse and even the expected length of a typical completed ray has 
not been found. There is no published paper by Gilbert on the topic; notes he supplied appear 
in a book by Noble, with due acknowledgement to Gilbert. Citations have typically attributed 
the notes to Gilbert (as we do, see [3]). 

A version of the model where the directions of growth were confined to two orthogonal direc- 
tions, vertical (V) and horizontal (H), was discussed by Mackissack and Miles [3]. A tessellation 
of the plane by rectangles results in their model. This structure too has not yielded to analysis, 
although when seeds are equally likely to be V or H the authors did provide an analytic ap- 
proximation (based on ideas of Gilbert) to the expected ray length, namely y/2/X. The merits 
of this approximation have not been evaluated in the literature to date. 
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The current paper arises from work done in 1997 by the second and third authors (Cowan 
and Ma). They obtained some analytic results for an even simpler V&iH-model, whereby the 
growth of eastward-growing rays is halted only by southward-growing rays (and vice versa). 
Westward and northward have the same reciprocity. A realisation of their tessellation is given 
in Figure [TJ 
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FIGURE 1. The Cowan-Ma (or half-Gilbert) rectangular tessellation when V-type and H— type 
seeds have equal intensities. 

Cowan and Ma placed a recurrence relationship (see ([I]) below) from their work on the internet 
[2J, though without proof. The background to this recurrence is as follows. 

Consider the isosceles right-angle triangle POQ in Figure E|a). Here \OP\ = \PQ\ = I- 
Suppose n seeds lie inside the triangle, uniformly and independently distributed; the figure uses 
n = 6. East or south growth of the rays is shown. Because of the blocking rules, only some of 
the rays reach the boundary of the triangle POQ. 

Cowan and Ma investigated the probability h n that no rays hit the boundary within the 
segment OP. This can also be interpreted as the probability that L, the final length of a test 
ray commencing eastward growth from O, is > £. 




Figure 2. Dia grams to assist the proof of the Cowan-Ma recurrence. 
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Their recurrence relationship for h n was as follows. For n > 1, 

n\ q ^ n s^ u 2 n ~ u ~ v h u h v (n - 1 + u - v)\ (n - 1 - u + v)l 
W n ~j2n)\^ ^ u\ v\(n-l-u-v)\ ' 

with ho = 1. Here q is the proportion of seeds which grow horizontally. The recurrence does not 
involve £, so h n does not depend on t — as is obvious from the scale invariance of the problem 
posed by Figure [2^a). 

This recurrence is a useful analytic step, providing precise information on E(L) and F{£) : = 
Pr{L < £}. 

n>0 

from which we deduce (in an extended notation which includes A) that F\(£) = F\{\/~\£). Also 

roc 

E(L) = / [1 - F(£)]d£ 
Jo 

l, POO 

y[hL / (A^ 2 /2) n exp(-A£ 2 /2)^ 

/i n r(n + i) 



n>0 



— r 

ATT 



v Z/V n>0 

In Section [21 we report the proof used to derive the recurrence relationship §1$ and plot the 
probability density function of the random variable L. The plot has an extraordinary property, 
discovered when certain simulations of the full rectangular-Gilbert model done by our first 
author, Burridge PQ, were also plotted. The probability density function of the Cowan-Ma 
model with A = 2 was indiscernible from that of the full rectangular-Gilbert model with A = 1. 

Section [3] presents Burridge's simulation study, that has a very high level of accuracy, and 
discusses this surprise coincidence — which raises somewhat the profile of the Cowan-Ma model. 
As well as having interest in its own right as a tessellation model with tractable mathematics, 
the model provides approximations for the full-Gilbert rectangular model. For example, the 
Cowan-Ma model — which we also called the half-Gilbert model because it has half of the 
blocking mechanisms — provides a much better approximation for E(L) in the full model when 
compared with the Mackissack/Miles approximation, which is E(L) ~ a/2/A when q = \- 

In Section [J] our work pushes further the tractability of the half-Gilbert model finding; most 
notably, we find that the mean ray length when q = ^ is given by the formula: 

E(L) 



VX(r(!)) 2 

In both our simulation and analytical work we have employed Zuyev's concept of stopping set 
sequences [5] and the distributional results for the areas of these sets. To achieve the analytic 
results, we have incorporated a new concept into the analysis, the idea of dead zones which 
influence the formation of the next stopping set in the sequence. Our most complete analysis is 
for the balanced case, q= \, because some results become rather complicated when q^ \. The 
expected ray length in the latter case is reported, without proof, in the appendix. 

2. The Cowan-Ma recurrence relation 

We now prove ([I]) for general q. Obviously, Hq = 1 and h\ = \. When n > 1, we label the seed 
closest to OP as A. See Figure 2(b). If the distance from A to OP is denoted by the random 
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variable Y, it is easily shown that Y has probability density function gy{y) = 2n(l — y) 2n ~ 1 , < 
y < 1. Furthermore, the conditional probability density function of X := \AB\ given Y is 

. . . 1 
9(x\y) = , 0<x<l-y. 

i - y 

Denote the event that no rays hit OP by £ n . Then 

Pr{£ n |a;,y} = Prjseed A grows eastward and reaches B 

and no ray grows across the segment EA} 

= q Pr{no ray grows across AB and no ray grows across EA} . 

To evaluate the right-hand side, we partition the domain above EB into the three zones that 
are shown in Figure 2(b). We then consider the trinomial distribution by which the remaining 
(n— 1) seeds are allocated to these zones: u to ABC, v to EAD and the remaining (n — l — u — v) 
to ACQD. This leads, for each (u, v), to a rather pleasing representation of the problem into two 
problems self-similar to the original one. Continuing, using || • || as area, we write Pr{£ n \x,y} 
as 

" (n - 1)\ \\ABC\\ U \\EAD\\ V WACQD]^- 1 -' 11 



it! v\ (n- l-u-v)\ WEBQW"- 1 



x 



Pr{no ray grows across AB and no ray grows across EA\u, v} 

~ q hh u\v\{n-l-u- v)\ ((i^)-i 

Pr{no ray grows across AB\u} Pr{no ray grows across 

^ ln ^" (n - 1)! x 2u (l -x- y) 2v [2x{l -x- y )]n~i-u~v 

~ q ^ n ^ u\ v\(n-l-U- V )\ (l-y)2(n-l) K K ' 

u=0 v=0 v J \ tfj 

Unconditional on x and y, and with n > 1, 
/i n = Pr{£ n } 

Pr{£|x,y}#y(y)sr(a;|y) etedy 
(n - 1)! /i M ^ 



o Jo 

n— 1 n— 1— m 

y y 

^— ' ^— ' u! vl (n — 1 — u — v)\ 

u=0 n=0 

"1 /•!-» 



2n / / x 2u (l - x - y) 2t '[2x(l - x - y)]^ 1 ^^ dxdy 
Jo Jo 



'o Jo 

n-ln-l-u on _ 1 _ u _ JJ 



u! vl (n — l — u — v)! 

u=o u=o v ; 

2n /" / " x n - 1+u - v {l-y-x) n - l - u+v dxdy 
Jo Jo 

n— 1 n —l— u nn—u—v u u 

n\ q } > — : — — — x 

^— ' ^— ' u! u! (n — 1 — u — 

u=0 u=0 v ; 

fl 

(1 - y) 2n ~ 1 B(n + u-v,n-u + v) dy 

o 

n\ q n y^ u 2"-"-" /i M K (n - 1 + tt - tpl (n - 1 - u + p)! 
(2n)\ ^ u\ v\ (n — l — u — v)\ 

^ 1 u =0 v=0 v ' 
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FIGURE 3. The three solid curves are the probability density functions / for the final length of a 
typical eastward-growing ray in the Cowan-Ma model. Each is based on A = 2 and with three q values: 
j, i and We show later in the paper that: (a) each of these solid curves actually comprises two 
curves overlaid, the second being the curve from the full-Gilbert model, with A = 1; (b) the dashed lines 
are the probability density functions from Gilbert's heuristic 'mean field' analysis, also valid for both 

models. 



We augment this recurrence with the result ho = 1. This completes the proof of dU). We note 
that the sequence ho, hi, hi, ... commences 1, |, 3, ^jj, ... when q = \. 

The recurrence together with ([2]) can be used to plot f(£) := F'(t) against I for various values 
of q (see Figure [3]) . 

3. Simulation of the full rectangular-Gilbert tessellation 

Finding coefficients analogous to h n for the full rectangular model is a formidable task because 
of the complexity of the blocking effects. Lacking self-similar zones akin to those discovered by 
Cowan and Ma in their model, we have devised an efficient way of accounting for these effects 
by simulation. 

The analogue of the isosceles triangle used in Figure [2] is a square, rotated so that its diagonal 
AC lies east-west, as illustrated in Figure [U To study the growth of horizontal rays, we consider 
an iJ— type test seed located at the western corner of the square, marked A in the figure, and 
define: 

h n = Prjray from test seed A reaches B \ n seeds in the square}. 

The only seeds that can block the test ray lie in the western side of the square, but whether or 
not they do so depends also on the configuration of seeds in the eastern side. Seeds outside the 
square have no influence. 

By analogy with equation (j2j) the ray length distribution for the rectangular Gilbert tessella- 
tion is: 

F W = 1 -yjh„' 2Ap >"7'- 2A<2 », 

n>0 

from which F\(t) = Fi(v2A^) is deduced. 

An obvious method: The naive approach to estimating h n would be to repeatedly populate 
the large square in Figure 0] with n seeds (each independently of ii-type with probability q) and, 
each time, determine if the line AB is intersected. This can be accomplished using the following 
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recursive algorithm which decides if a ray, extending in compass-direction u G {— >, f , J,} from 
one seed s* will be blocked within a distance d The algorithm, block, outputs a logical value: 

, / * , x f true if ray is blocked 
block s , a, it) = <„_, ■ , t i i i 

v ' I false if ray is not blocked. 

Whether or not s* is blocked within a distance d depends only on the configuration of seeds 
within a square of diagonal 2d along which its produced ray travels. Let the compass-direction 
of this ray be u, and let us denote by A(s*,d,u) the isosceles triangle which forms the half of 
the square closest to s*. Let the type (H or V) of seed s be t(s). The algorithm block(s*, d, u) 
runs as follows: 

for all s e A(s*,d,u) do 
if t(s) ^ t(s*) then 

compute the perpendicular distance, d s , and compass-direction, u s from s to s*'s pro- 
duced ray. 

if block(s, d s , u s ) = false then 

return true 
end if 
end if 
end for 
return false 




FIGURE 4. Only particles within the large square (whose diagonal is AC) can influence the event 
that an H— type ray starting at A does not reach B (due to intersection of the line-segment AB by 
vertical rays). The role of the smaller shaded square is described in the text. 



For example, if s* is the H-type seed at A, then the computer programme calls block(s*, £, — >). 
This invokes recursive calls to block for every V-type seed in the left isosceles triangle (until 
a true value is returned by the call). In Figure |U the shaded region with a F-type seed s at 
the top shows a square that is investigated by one of the recursive calls, specifically by the call 
block(s, d s , 4), where 2d s is the diagonal length of the shaded square. 

In principle, we can conduct this simulation for each n up to (say) 300. For each n, we would 
generate the seeds in the square (with diagonal AC) N times, where N would be very large. An 
estimate of h n ,, < n < 300, is thereby generated for Zf-type rays. Then, if q ^ ^, we would 
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repeat the whole procedure for V-type rays. It is a lengthy process, despite the potential saving 
if an early-tested seed s returns true — implying that others don't have to be tested. 




(a) (b) (c) 

FIGURE 5. Nested stopping sets are shown in (a) and (b). In (c), the 'efficient algorithm' is at step 3 
and at least one more step will be needed before we see a ray in step n crossing the dashed half-diagonal 

of S n +1. 



Stopping sets: To shorten the task, we have devised a method based on stopping sets (a 
concept defined by Zuyev [5] and amplified in [6]). Consider the unbounded quadrant that lies 
between the half-lines y = x and y = —x, with x > 0, partly shown in Figure [5]^a). A stationary 
Poisson process of seeds with intensity A exists in the quadrant. A triangular set whose eastern 
boundary is vertical and western vertex is the quadrant's apex is gradually expanded, stopping 
briefly whenever its boundary hits a seed — before continuing its expansion. The set stopped 
by the kth seed encountered is called S^. This process creates a nested sequence of random sets. 
We denote the area of S% by E\ and the areas of the region \ Sfc-i by Ek, k > 1. Another 
nesting arrangement is shown in Figure EJb), this time with squares and a different ordering of 
the seeds. 

Zuyev showed, among other things of a more general nature, that any expanding domain 
constructing a nest of compact sets in the manner described above — through a sequence of 
stops caused by seed hits — creates areas E\ , E2 , E3 , . . . which are independent and distributed 
exponentially with parameter A. The domain might have a complicated geometry because 
the expansion rule is allowed to depend on the seeds that it contains (and, being closed, this 
includes seeds on the domain's boundary). In the two examples of Figure O the expansion rule 
is straightforward and doesn't depend on the internal seeds. 

Most importantly for the validity of Zuyev's distributional results, neither the expansion 
rule nor the stopping rule for Si should depend on seeds outside the expanding domain. This 
prohibition plays two roles. 

• it helps establish that E\ is exponentially distributed; 

• it also allows one to say that the point process of seeds outside the stopping set S\ is 
still a stationary Poisson point process with unchanged intensity given the information 
within S\ (a notion formalised by Theorem 2 of [6]). 

This allows the argument to be extended sequentially to E2, E3, ... and S2, S3, 
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We also note that Zuyev's results are not guaranteed if randomisations apart from the Poisson 
process of seeds affect the growth and stopping. No such complication occurs in this section of 
our paper, although we must address the issue later in Section |H 

Stopping sets constructed in this way have other properties. The ith seed s, is uniformly dis- 
tributed on the growth frontier of Si and the seeds s\, S2, s n are uniformly and independently 
distributed in the set SVi+i- Those of Figure 0(b) have a property that no other nesting has: 
if Si is V— type, then whether or not it reaches the east-west diagonal depends only on seeds 

Sl,S 2 , ...,Si-l. 

Efficient algorithm: In the context of Figure 0(b) with its nesting of squares, the latter 
property says that the ray growth just within S n +\ from the seeds s±, S2, ■ s n provides a sample 
of the problem that interests us — giving a true or false datum on whether a test ray is blocked 
before it traverses across half the diagonal of S n +\. (See the illustration for n = 3 in Figure 
Efc).) This datum contributes to the estimation of h n . Importantly, as we show below, if the 
datum is true, then we can add a true datum for the estimation of all hj,j > n — without 
further computational effort. 

We start with the unbounded quadrant empty of seeds, then place an if-type test seed at 
the apex of the quadrant. We generate the exponentially distributed areas E\ and E2 and so 
construct the squares S\ and 52 expanding from the apex. We randomly select (uniformly) a 
seed point s\ on the growth frontier (eastern sides) of the inner square, S\. Because of the 
properties discussed above, this is equivalent to choosing the point uniformly within the outer 
square S2. If this seed grows a vertical ray that intersects the diagonal, let the distance of the 
intersection point from the apex be X±. If not, set X± = 00. 

Let Ai denote \\Si\\, the area of Si, and £ n denote the event that the line from the test 
seed reaches the centre of a square populated with n uniformly distributed seeds. Obviously 
h n = Pr{£ n }. 

If X\ < {A2/2) 1 / 2 then the simulation ends. There is no need to generate more nested squares 
in order to simulate the events £ n , n > 1 because we know that the half diagonal of every 
subsequent square will be crossed at X\ < {An/2) 1 / 2 . Seeds on the boundaries of subsequent 
squares cannot influence this. If the first seed does not cross the diagonal, or crosses such that 
X\ > {A2/2) 1 / 2 , then we draw S3 and pick a point S2 on the boundary of the second nested 
square S2. We check if S2 intersects the diagonal, accounting for any possible blocking effects 
from si by using the algorithm block. If so, we let the distance from the apex to the closest 
intersection point be X2, which will be < X±. If X2 < (A3/2) 1 / 2 then the simulation ends. If 
not, we add another square S4 and seed S3 — reaching the situation in Figure [5](c) — and so 
on. We keep repeating the process — adding another seed and using block on that seed — until 
block indicates that the latest half-diagonal has been hit. We then record that the event £ n 
fails to occur for this and all higher values of n. The entity h n for eastward growing rays is 
the fraction of times that £ n occurs over many simulations. If q 7^ ^, the complete protocol is 
repeated with q replaced by (1 — q) to give results for southward growing rays. 

To estimate the h n , N = 10 9 simulations were performed, requiring a running time of approx- 
imately one hour on a modern PC. When q = \, the largest number of nested squares created 
before the simulation terminated was 917, which occurred once, and the second largest number 
was 727, which also occurred once. The mean number of squares created before termination was 
5.25. In the (q = |) case, the estimate of expected length of each line produced from a seed 
was: 

(3) E(L) = 1.467535 (0.000029) 

where the bracketed number is the standard error, calculated with due regard to the positive 
covariance between our estimators of h n and h n+ jt, k > 0. 

Remark 1 : Our accurate estimate of the h n values allows the probability density function of 
the ray length to be calculated. Because the two ray lengths coming from a particular seed are 
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independent, the standard convolution method leads to an estimated distribution of the total line 
length arising from a typical seed. Mackisack and Miles [4j claim that these two ray-lengths are 
not independent, but we disagree. 

The coincidence: We found a remarkable similarity between the probability density func- 
tions in the half rectangular Gilbert model and the full rectangular Gilbert model when the 
intensity of seeds in the former case was twice that of the latter case. Indeed the plots were 
almost indistinguishable, so Figure [3] effectively displays both / and f for various q, with A = 2 
or A = 1 respectively. 

We are mindful that the simulated results have sampling error, albeit small. So we asked the 
question: are the two distributions F and F mathematically equal — or just approximately so? 
To answer this in the (q = i) case, we performed some rather tedious exact calculations (details 
omitted) which yielded: 

ho = 1; hi = |; h 2 = jr;', h 3 = ^. 

We then expanded both F and F as Taylor series about the origin. 

F{£) = h + (h x - h )£ 2 + l(h + h 2 - 2/nK 4 

+ i(3/ii - 3/i 2 - ho + h 3 )£ 6 + o(£ 7 ) 

= l-\e + \^-§-^ + {f). 

F(£) =h + 2(hi - h )£ 2 + 2(h - 4hi + h 2 )£ 4 

+ |(3hi - 3h 2 - h + h 3 )£ 6 + o{£ 7 ) 

= l-\e + \£±-^ + o(£ 7 ). 

We see that these exact series differ slightly in the fourth term, so F and F are not mathematical 
equal. 

"Mean field" analysis when q = |: Gilbert's original "mean field" analysis, which 
was adapted by Mackisack and Miles [I] to the (q = ^) rectangular case, involved the rough 
approximation that ray ends (there being two per seed) were uniformly spread across the plane. 
With this assumption, it was possible to approximate at time t the expected number of ray ends 
lying within a small distance Sx of rays that would block the growth of these ends within the 
next 5t. 

Mackisack and Miles analyzed the {q = \) full model using two quantities; R(t), the expected 
total length of rays per unit area; G(t), the expected number of growing ends per unit area. 
Recounting their work, these quantities are related exactly by R = G, assuming unit growth rate, 
and heuristically in the full rectangular case by G ~ -^RG, with initial conditions R(0) = 

and G(0) = 2A. Solving these differential equations, they found that G(t) ~ 2Asech 2 y^f t. If L 
is the final length of a test ray in their full Gilbert model, then: 

(4) Pl{L>e)= m = m^f lt 

The expected L when q = | is therefore approximated by y / 2/A = 1.41421 at A = 1. This is not 
especially close to the value shown in ([3]). The solution for R was R(t) ~ 2\/2A tanh(y / 2~7A i), t > 
0. 

We have modified the analysis in [3] to deal with the {q = \) Cowan-Ma model. We put 
G ~ — \RG since each of the four directions of growing lines can only be blocked by one other 
line type. Solving the new equation pair, we find that the number of growing lines per unit area 

at time t for the half model is : G(t) w 2Asech 2 y / ^. Also R(t) m 2 V / 2A tanh(y / 27A t), t > 0. 
Furthermore (|4|) still holds, with G replacing G. So, setting A = 2 in the half system and A = 1 
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in the full system we obtain identical approximations to the probability density function for ray 
length: 

f 2 (£) « V2 sech 2 -^= tanh -= » fi(^). 
v2 v 2 

The expected ray length is: E(L) ~ y/2. 

So we have shown that the mean field approximations in the two models are equal, when 
q = Indeed, our analysis for q ^ i, developed in the next sub-section, shows that the two 
approximations are also equal when q ^ ^. 

Mean field analysis when g / ^: When the intensities of H- and V-type seeds are 
not equal, the rays of east-growing and south-growing have different length distributions. So 
a system of four differential equations and four initial values is needed, in variables (for the 
half-Gilbert model) R± and .R_>. 

R+(t) = G^(t) G^(t) « -R±(t)G->(t) 
R x (t) = G x (t) G±(t) » -R^(t)G x (t), 

combined with: 

iU(0) = Ri(0) = 0; G_>(0)=gA; G+(0) = (1 - q)X. 

Replacing ~ with = and eliminating Gh and Gy, the differential equations become 

= -R±(t)R^(t) 
R^t) = —R-+(t)R±(t), 

augmented by 

iU(0) = R^O) = R^(0) = qX R^O) = (1 - q)X. 

We have only been able to solve this coupled system in series form and, even then, with no 
general term recognised. Using the abbreviations Q := qX and P := (1 — q)X, 

p M ®+ PQ + , PQ(3P + Q) 5 PQ(15P 2 + 16PQ + 3Q 2 ) ^ 7 

^ () = n + — si — 7i 

PQ(105P 3 + 241P 2 Q + 135PQ 2 + 15Q 3 ) 9 
[p) + gj * - 

with G_>.(i) being R->(t) (easily calculated from [5]). A Mathematica routine to compute as many 
terms as required is available from the authors. For R^ and G^, simply interchange P and Q. 
Note that west-growing rays have results identical to east-growing — likewise north and south 
results are identical. 

For the full-Gilbert model, the equations are very similar, but cast in terms of the four variates 
Gy,G}j,m,y and R//. 

R H (t) = G H (t) G H (t) » -K v (t)G H (t) 

Ry(t) = Gy(t) Gy(t) « -R H (t)Gy(t), 

combined with 

Rh(0) = Ry(0) = Gjy(0) = 2qX Gy(0) = 2(1 - q)X. 

This leads to a solution for R#(t) equal to the right-hand side of ([5]), but with Q = 2qX 
and P = 2(1 — q)X. Thus it becomes obvious that R#(i) with A = 1 equals R^(t) with 
A = 2. Likewise for the other linked pairs of variables! Therefore, when q ^ \, the two ray 
length distributions (for H and V rays) for the full model having intensity A are equal to the 
corresponding ray length distributions for the half-Gilbert model with seed-intensity 2A. All of 
these entities are, of course, only approximate solutions to the true Gilbert models. 
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Figure [3] shows that their value as approximations for the full-Gilbert model is quite good, 
but not nearly as good as the analytic answers adopted from the half-Gilbert model. In the 
last section of the paper, we provide more of these answers, demonstrating that the half-Gilbert 
model of Cowan and Ma is encouragingly tractable. 

4. Stopping sets and dead zones in the half-Gilbert model. 

It is possible to use the stopping-set concept to find exact expressions for the first, second 
and in principle higher moments of the ray length in the Cowan-Ma model. The balanced case, 
q = ^, is easier to describe — and that is now our focus. We give some results for the general 
case in an Appendix. 

A different construction of stopping sets: Suppose that a stationary Poisson process of 
intensity A exists in the plane, with seeds marked either H (east growing) or V (south growing) 
with equal probability. In Section [3] we have described how a nest of Zuyev's stopping sets is 
created when the growth frontier of an expanding domain hits the seeds. For the Cowan-Ma 
model, the seeds that are relevant for an east-growing test ray commencing at O in Figure [2{a) 
is the shaded region in that figure — or, more precisely, the unbounded octant lying between 
y = x and y = 0, with x > 0: we call this region, the initial live zone. 

As before we start by expanding a domain — an isosceles right angle triangle in this case (see 
Figured]) — into the live zone, stopping when it hits the first seed s\ whose coordinates relative 
to O are (x\,yi). This creates a domain Si with area E\ that is exponentially distributed. If 
s\ is V— type, then it will provide the ray that blocks the test seed; thus L = x\ and no other 
seeds need be considered. 

Alternatively if s\ is H— type, then, instead of growing Si (retaining its shape as an isosceles 
right-angle triangle and constructing the familiar Zuyev nest of stopping sets), we introduce 
a significant modification. We remove a part of the live zone: a 'dead zone' labelled D\ (see 
Figure [6]) which has now become irrelevant, as we shall soon see. 

As S\ U D\ has been constructed without drawing upon any information taken from outside 
Si U Di, the point process in the remaining region (the new live zone) is still a Poisson process 
with unchanged intensity given the information within Si U Di — as explained in Section [3j 

We now grow a trapezium whose left-hand side located at x = xi has length y = y\. The 
trapezium expands until its right-hand side first hits a seed S2 (in the new live zone). The 
stopping set formed is called S2. It has an exponentially distributed area Ei- 

We proceed in this way, forming a sequence of stopping sets (illustrated in Figure [6]) which, 
unlike those in Section[3j do not form a nest. They do, however, have independent exponentially- 
distributed areas and are part of a recursive structure which we can exploit. It is also important 
to note that the first V-type seed will provide the ray which blocks the test ray. Without our 
introduction of dead zones, a complicated algorithm rather like block would be required to 
check if a y-type ray actually reaches the path of the test ray. 

Remark 2: Why is it that no seed in Di can influence the distance L travelled by the test 
seed; dead zone V -types will either be blocked by an east-growing ray within the live zone or, 
if they are not blocked, the test ray must have been intersected at an earlier point. Dead zone 
H-types can never be in a geometrical position to block live zone V -types. The same line of 
reasoning applies to dead zones D2,D%,.... 

Remark 3: In Section 3 we mentioned that extraneous randomisations, those not solely 
dependent on the Poisson point process, might invalidate the key results from the stopping-set 
theory. There is no such problem here with the stopping set Sk itself, but we note that D^, k > 1, 
depends for its existence on an extraneous random feature — namely the H or V mark of seed 
Sf-. This does not invalidate our comment above that the point process in the current live zone 
outside Sk U is unaffected by the information in Sk U D^. For one thing the seed marks are 
independent of each other and of the point process. Furthermore, we only stop constructing dead 
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zones when we have no further need to observe the process at all. So the extraneous random 
feature is not operative in our analysis. 

The recursive structure commencing with a generic live zone: Suppose that we begin 
observing the process when the live zone has left boundary of height y and when we are about 
to construct S n . In Figure El we draw the case n = 2. The probability density function for 
the length, r, of S n 's base, conditional on the height y of its left boundary, follows from the 
exponential distribution of 5 n 's area E n . It is therefore: 

f(r | y) = A(r +y)e-^ r2+2ry \ 

If the stopping seed s n for set S n is V-type, then its south ray will be the first to intersect the 
test ray and the process ends. Otherwise, another dead zone is created and further trapezoidal 
stopping sets are formed until a V-type is met. 

Let X be the random variable equal to the horizontal distance covered by stopping sets until 
the process comes to an end. The density function of X, conditional on y will be: 



(6) g(x\y) = ^ 



(x + y) e -^ x2 + 2x y) + J™ e -U r2 + 2r v) g( x - r\u)du^j dr 



where g{x\y) = if x < 0. The first term in the square brackets accounts for the case where 
the first seed is V-type, and the second term for the case where it is H-type and the process is 
effectively re-started with a different boundary condition having already covered some horizontal 
distance. We have taken q = |, but the analysis can be carried out for general q, producing a 
more complicated result. Note that the ray length probability density function is g(x \ 0). 
We define the moments of the conditional density: 

/•oo 

IJ-n(y) = / x n g{x\y)dx. 
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As mentioned before we will here compute E(L) = /xi(0) and E(L 2 ) = /i2(0), which from equation 
([!]), satisfy: 



(7) 



(8) 



e 2" 2 fii(u)du. 



Our strategy is to find fJ.i(y) and ^{y) up to an arbitrary constant, and then to determine the 
constant using equations ([7J) and ([8]). The first part of this process is most easily achieved by 
making use of the moment generating function: 



M t (y) 



which, from equation (0) satisfies 
M t (y) 



1 1 / 7T (Ay~t) 2 

— I — \ / — e 2A { erfc 

2 2 V 2A 



Ay-* 



t + A / M t {u)du 



+A 



erfc 



Xu — t 



Mt(u)du 



This integral equation may be reduced to the differential equation: 

(Xy-t)- 



d 2 M t „ , N dAf t A_ A 



ciy 2 " ; 2 M * 2' 

Expressing the left hand side as a series in t, and collecting coefficients of t and t 2 we obtain 
differential equations satisfied by n±(y) and ^{y) '■ 



(9) 
(10) 



fi'l(y) - Xy^[(y) - ^/ii(y) = 



Clearly we must solve for fii(y) first. 

The first conditional moment: Making the change of variable z = (^)^y in equation 
we obtain: 

d 2 fii dfij 



dz 2 



2z-p -ui = 0. 

dz 



If the coefficient of \i\ were a positive multiple of two, this would be Hermite's equation, solved 
by Hermite polynomials. Since this is not the case, we seek a series solution [7]: 

oo 

v\{y{ z )) = J2 anZn 



and obtain the recurrence relation: 



«n+2 



n=0 



2n+l 



(n + l)(n + 2) 



This leads to the general solution: 



Hi(y{z)) =a Ma±,z 2 ) + a lZ M(f,f,z 2 ), 



where M is Rummer's Function 



M(a,b,z) = 



n=0 



{a)nz n 
(b) n n\ ■ 
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Here we have used the Pochhammer symbol, defined by: 

(a) n = a(a + l)(a + 2). ..(a + n — 1), (a) = 1. 

The Kummer's functions diverge as z — > oo, but we know that fii(y(z)) — > in that limit. This 
apparent paradox is resolved by noting that the two independent parts of the solution may be 
combined to form a Kummer's function of the second kind [8], defined by: 



7T 

U(a, b, z) 



sin7r6 



M(a,b,z) ^_ b M(l + a-b,2-b,z) 



T(l + a-b)T(b) T(a)T(2-b) 



which tends to zero as z — > oo. In terms of this function, the general solution is: 

Hx{y(z)) =Ax M(i, i, z 2 ) + B x U{\, \,z 2 ). 

It must be the case that A = in order to capture the right asymptotic behavior so, restoring 
the original variable y, the conditional moment must have the form: 

(11) n 1 (y) = BxU(\,±,±y 2 ). 

It now remains to compute B. We do this by substituting (llip into equation ([7|). Making use 
of the result: 

\ erfc(u)[/(i, \y)du = ^[r(i) - v^r(f )], 
Jo 7r 

together with r(|)r(|) = V^ir we find that: 



B J7T 



We have now found fii(y), which gives us a compact analytic expression for the expected ray 
length: 



E(L) = m (0) 



7T 



VXr(f) 

7T 

= VX(r(f)) 2 

2.0920992 

For comparison, using the first 200 coefficients from Cowan and Ma's recurrence we obtain 
E(L) ~ 2.0920987 when A = 1. As we discovered earlier, when A = 2, the half model provides an 
approximation to the full model, having similar but simplified blocking effects and identical mean 
field behaviour. For this choice of A we obtain the exact half model result E(L) = 1.479337560 to 
7 decimal places, which differs from the accurate full model result (1.467535) by 0.7%. Compared 
with the mean field prediction: E(L) ~ \/2, which differs from the full model by 3.6% this is a 
much closer approximation. 

The second conditional moment: As for the calculation of fii, we make the change of 
variable z = (^)iy, but this time in equation (|10p . obtaining 

dz 2 dz w "Ar(f) U ' 2 ' h 

where we have used the differential property [8] : U'(a, b, z) = —a U(a + 1,6 + 1, z). We know 
the homogeneous part of the general solution to (jlOp . so it remains to find a particular solution. 
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We do this using variation of parameters, and begin by making the definitions: 

f 1 (z) = M(l 1 1 ,z 2 ) 
f 2 {z) = U{\,\,z 2 ). 

The function M has the differential property [8]: M'(a, b, z) = | M(a + 1, b + 1, z) which allows 
us to compute the Wronskian: 

W(z) = h{z)f! 2 (z)-h(z)f[{z) 
z 
2 

We now define 

G(z,t) 



5 3 z 2. 



- [M(i, I, z 2 )U{\, §, z 2 ) + 2U{\, i, z 2 )M(l I 
f2(z)f 1 (t)-f 1 (z)f 2 (t) 



in terms of which the particular integral is: 

f I \ - ^ 

U[z) ~ Ar(i) 

Discarding the divergent part of the solution, and restoring y, we have that: 

J2tt 



W(t) 



G(z,t)tU(Ht 2 )dt. 



fl2 ( y ) = CxU(\,±,±y 2 ) 



Ar(|) J ziy) 



G(z(y),t)tU(Ht 2 )dt, 



where C is an as yet undetermined constant. We find it by substituting our expression for ^{y) 
into equation (|8|). Making use of the numerical integral: 



K 



erfc(z) 



o 



G(z,t)tU(Ht 2 )dt 



dz 



we find that: 



0.343146 



C 



ttK 



r(!)A\r(f) 



+ 2V2 . 



Noting also that / p (0) = j we have the final result that: 



E(L 2 ) = M2 (0) 
1 



ttK 



+ 2V2) [/(n,o) + / p (o) 



r(f)A Vr(|) 
7r 3 / 2 K + 2r(f) (v^F + r(|) 



3\ 3 



Ar(f) 



6.37688 
A 



For comparison, using the first 200 coefficients from the Cowan-Ma recurrence we obtain E(L) 
6.37686 when A = 1. 



5. Concluding comment 

Gilbert's tessellation is notoriously difficult to analyse, and even the rectangular version stud- 
ied by Mackisack and Miles remains entirely without analytical results. In this paper we have 
shown that the simplified rectangular model of Cowan and Ma, with only half of the blocking 
rules of the Mackisack and Miles model, has a number of tractable properties. As such, it is the 
only Gilbert-style model, we believe, which has yielded any analytic results. 
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Appendix: Expected length in the half model when q ^ ~ 

If q is the proportion of seeds growing horizontally in the half model, then equation ([6]) 
becomes: 

rr+y 



g(x\y) = (1 - q)\{x + y) e -^ x2+2xy) + qX e -^+2ry) 

Jo 



g{x — r\u)du 







dr. 



E(L H ) 
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FIGURE 7. Expected horizontal length in the half model as a function of q, the proportion of 
horizontal rays. The seed density is A = 1. 



The first moment of g(x\0) may be found by similar methods to those employed in the q 
case. The expected length of a horizontal ray is found to be: 

-l 



V2- 



(i - 1) Gli (i 



013+1 
1 -1 



2 9+ Wr(i 



where G is Meijer's G-Function [9]. Figure [7] illustrates the function E[L#] for the case A = 1. 
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